SIM202 : HOME SAVER

Cyril NEFZAOUI & Bertrand PATUREL

Présentation des codes et des résultats

Sommaire :

Visualisation des données

Imputation par continuité

Modèle linéaire

GAM

GBM

XGBoost

Imputation VS prévision pure

Permutation

Complétion partielle

Initialisation

On récupère les deux jeux de données train (jeu de données d'entrainement) et test (jeu de données de submission).

In [13]:
setwd("C:/Users/patur/OneDrive/Bureau/ENSTA/ENSTA Cours/Projet R/TP")

rm(list = objects())#comme clear all, met a jour des var

train=read.csv("Data/train.csv",header=TRUE) #on charge le tableau d'entrainement train
test=read.csv("Data/test.csv",header=TRUE) #on charge le tableu de résultat test

Visualisation des données

On regarde tout d'abord a quoi ressemble notre jeu de données Appliances en fonction du temps.

In [14]:
library(ggplot2)
ggplot(train,aes(Posan,Appliances))+geom_point()+xlab("Temps")+ylab("Consommation")
Warning message:
"package 'ggplot2' was built under R version 3.6.3"

Puis l'on s'interresse au jeu de donnée test. On souhaite comprendre quelle est la plage des relevés de cette table et la comparer aux dates du jeu de données train.

In [15]:
library(xts)

Date_train=as.POSIXct(strptime(train$date,tz="GMT","%Y-%m-%d %H:%M:%S")) #conversion de la date de train
train.xts=xts(train$Appliances,order.by = Date_train) #conversion de la liste train en série temporelle

end=end(train$Posan)[1]
set=which(test$Posan==train$Posan[end])
l=end(set)[1]
first=set[l]
last=end(test$date)[1]

Date_test=as.POSIXct(strptime(test$date,tz="GMT","%Y-%m-%d %H:%M:%S")) #conversion de la date à prévoir de test 
num=last-first
c=rep(200,num)
test.xts=xts(c,order.by = Date_test[(first+1):last]) #conversion de la liste à prévoir du test en série temporelle

res.xts=cbind(train.xts,test.xts)
plot(res.xts) #en noir ce que l'on connait en rouge ce que l'on doit prévoir
Warning message:
"package 'xts' was built under R version 3.6.2"Loading required package: zoo
Warning message:
"package 'zoo' was built under R version 3.6.3"
Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric


Attaching package: 'xts'

The following objects are masked from 'package:dplyr':

    first, last

Conclusion : On suppose donc qu'il y avait un jeu de données complet initialement qui à été répartie en deux jeux de données, l'un train et l'autre test

Avant tout traitement de train et test, nous avons rajouté ces data frame d'une colonne isintest, dont les coefficients valent 1 ou 0 selon que l'élèment associé appartient au train ou au test. On a ensuite fusionné train et test dans une large dataframe que l'on appelle total. On effectue un tri de total par ordre chronologique croissant, ce qui nous permet de visualiser la répartion de certaines valeurs du test dans le train.

In [16]:
#### On reconstitue le jeu de données initial 
library(mltools)
library(data.table)

#On créer une variable qui vaut 0 si le relevé est issue de train, 1 si c'est dans le test
yes_is_in_test=rep(1, times = length(test$date)) #On créer la varaible de 1
test$is_in_test=yes_is_in_test #On attache cette variable au test
not_in_test=rep(0, times = length(train$date)) #On créer la varaible de 0
train$is_in_test=not_in_test #On attache cette variable au train

#On reconstitue le tableau initial qui a donné forme aux deux jeu de données
train_dt=setDT(train, keep.rownames=FALSE, key=NULL, check.names=FALSE)
test_dt=setDT(test, keep.rownames=FALSE, key=NULL, check.names=FALSE)
total= merge(train_dt,test_dt,all=TRUE) #On recolle nos deux tableaux ensemble
#length(which(is.na(total$Appliances))) les valeurs manquantes des appliances du test à compléter figure comme des NA, il y en bien 5771 
#length(total$lights) vaut 19735

#Puis on réeordonne la série totale en fonction des dates
total$date <- as.Date(as.character(total$date),tz="GMT",format="%Y-%m-%d %H:%M:%S")
total=total[order(total$date,total$NSM, decreasing=FALSE),]
total=total[1:(19735-(5771-4654)),] #On prend seulement les dates à imputer
Warning message:
"package 'mltools' was built under R version 3.6.3"Warning message:
"package 'data.table' was built under R version 3.6.3"
Attaching package: 'data.table'

The following objects are masked from 'package:xts':

    first, last

The following objects are masked from 'package:dplyr':

    between, first, last

Traitement des données

On cherche en priorité à savoir ou sont les données manquantes dans les jeu de données.

In [1]:
#summary(train) Nous donne les colonnes avec les NA, on remarque donc qu'il y a seulement les colonnes RH_6 et Visibility
#summary(test) Nous donne les colonnes avec les NA, on remarque donc qu'il y a seulement les colonnes RH_6 et Visibility

On remarque que les colonnes RH 6 et Visibility qui correspondent à la visibilité et à l'humidité sur le facade extérieure de la maison ont de valeurs manquantes. Affichons ces deux variables pour comprendre pourquoi il y a des données manquantes :

In [49]:
ggplot(data=train) + aes(x=as.numeric(date),y=RH_6)+ geom_point() + xlab('relevés')+ylab('Visibilité')+ggtitle('Mesure de la visibilité en fonction des relevés')
Warning message:
"Removed 1541 rows containing missing values (geom_point)."

On comprend donc que les valeurs manquantes sont dus à la saturation des capteurs de l'humidité et de la visibilité lorsqu'il y a des conditions extrêmes.

In [50]:
index=which(is.na(train$RH_6))
#mean(train$RH_6[index-1],na.rm=T) vaut 12 
#var(train$RH_6[index-1], na.rm=T) variance faible 
summary(train$RH_6[index-1])


index=which(is.na(train$Visibility))
#mean(train$Visibility[index-1],na.rm=T) proche de 10 
#summary(train$Visibility[index-1])
print(train$Visibility[index+1],na.rm=T,na.print="Saturation")
#train$Visibility[index+1]
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  10.00   10.43   11.26   12.29   13.51   19.20    1494 
 [1] Saturation Saturation Saturation Saturation Saturation Saturation
 [7] Saturation Saturation Saturation Saturation Saturation Saturation
[13] Saturation Saturation Saturation Saturation Saturation Saturation
[19] Saturation Saturation Saturation Saturation Saturation Saturation
[25] Saturation Saturation   10.66667 Saturation Saturation   19.00000
[31] Saturation Saturation   11.50000 Saturation Saturation Saturation
[37] Saturation Saturation   10.33333 Saturation Saturation Saturation
[43] Saturation Saturation Saturation Saturation Saturation Saturation
[49] Saturation Saturation Saturation Saturation Saturation   10.00000
[55]   22.66667 Saturation Saturation   16.66667 Saturation Saturation
[61] Saturation Saturation Saturation Saturation Saturation Saturation
[67] Saturation Saturation Saturation Saturation Saturation Saturation
[73] Saturation Saturation Saturation Saturation   11.83333 Saturation
[79] Saturation Saturation Saturation Saturation   28.00000   11.33333

On choisit alors deux méthodes pour compléter les valeurs saturées.

  • Première méthode : Pour _RH6 on fixe simplement à 10 les valeurs manquantes
  • Deuxième méthode : Pour Visibility on prend la moyenne pondérée aux bords non saturés : VtNA=kj+kVtj+jj+kVt+k(Vt) représente la série des visibilités.
In [17]:
sequence1=which(is.na(train$RH_6))#Nettoyage des données de RH 
for (i in sequence1)
{train$RH_6[i]=10} #Pour les NA on prend comme valeur 10

sequence1=which(is.na(test$RH_6))#Nettoyage des données de RH 
for (i in sequence1)
{test$RH_6[i]=10} #Pour les NA on prend comme valeur 10


sequence2=which(is.na(train$Visibility))
for (i in sequence2)
{
  k=0
  j=0
  while (is.na(train$Visibility[i+k])==is.na(train$Visibility[i]))
  {k=k+1}
  while (is.na(train$Visibility[i-j])==is.na(train$Visibility[i]))
  {j=j+1}
  train$Visibility[i]=0.5*(train$Visibility[i-j]+train$Visibility[i+k])
}#On prend comme valeur pour les NA la moyenne des bords 

sequence2=which(is.na(test$Visibility))
for (i in sequence2)
{
  k=0
  j=0
  while (is.na(test$Visibility[i+k])==is.na(test$Visibility[i]))
  {k=k+1}
  while (is.na(test$Visibility[i-j])==is.na(test$Visibility[i]))
  {j=j+1}
  test$Visibility[i]=(k/(j+k))*(test$Visibility[i-j])+(j/(j+k))*(test$Visibility[i+k])
}#On prend comme valeur pour les NA la moyenne des bords
Error in (j/(j + k))(test$Visibility[i + k]): tentative d'appliquer un objet qui n'est pas une fonction
Traceback:

Imputation statistique des valeurs maquantes par continuité

Nous venons de mettre en évidence qu'il y avait deux parties, une partie comprenant l'imputation statistique et l'autre de la prédiction. Dans cette partie nous nous concentrons sur une méthode concise.

Notons (At) la série de toutes les Appliances, (Ate) celle du test à compléter et (Atr) celle du train. Pour chaque valeur Appliances (Ate) manquante dans le test, on regarde les deux valeurs de train At+kr et Atjr encadrant et on effectue une moyenne pondérée de ces deux valeurs. Puis on améliore cette aproximation en attribuant un poids différent à ces deux Appliances encadrantes. Ainsi plus la valeur manquante à compléter du test est proche temporellement de son bord gauche (resp droit) plus le poids attribuer au bord gauche (resp droit) est important.

On utilise donc la relation suivante : Ate=kj+kAtjr+jk+jAt+kr

In [18]:
library(ggplot2)
res=test$rv1
p=1 #Initialisation
sequence3=which(is.na(total$Appliances))
for (i in sequence3)
{
  k=1
  j=1
  while (is.na(total$Appliances[i+k])==is.na(total$Appliances[i]))
  {k=k+1}
  while (is.na(total$Appliances[i-j])==is.na(total$Appliances[i]))
  {j=j+1}
  total$Appliances[i]=(k/(k+j))*(total$Appliances[i-j])+(j/(j+k))*(total$Appliances[i+k])
  res[p]=(total$Appliances[i])
  p=p+1
}

ggplot(data = total)+geom_point(size=0.01)+aes(x = date,y = Appliances)

Le modèle linéaire

Choix du modèle linéaire

Ici on teste tout d'abord plusieurs manière d'effectuer un STEPAIC. On choisit finalement la méthode forward qui nous laisse plus aisément le choix des variables.

In [19]:
library(MASS)
level=13000 #On choisit comme donnée d'entrainement les 13000 premières valeurs du train

# On décide de retirer les variables redondantes et non numérique qui n'auraient pas de sens dans une regression
Data=train[1:level,-1] #On enlève la date qui ne peut être lu
Data=Data[,-33]
Data=Data[,-30]
Data=Data[,-30]
Data=Data[,-30]
fitstart=lm(Appliances~1,data=Data)
fitall=lm(Appliances~.,data=Data)

#Trois manière de trouver le modèle linéaire minimisant l'AIC
step.modelf=stepAIC(fitstart,direction="forward",scope=formula(fitall),trace=FALSE)
step.modelb=stepAIC(fitall, direction = "backward",trace=FALSE,steps=50)
step.model=stepAIC(fitall, direction = "both",trace=FALSE,steps=50)

#step.modelf$coefficients
#step.modelb$coefficients
#step.model$coefficients

# On décide de retirer les variables redondantes et non numérique qui n'auraient pas de sens dans une regression
Databis=train[level:13964,-1] #On enlève la date qui ne peut être lu
Databis=Databis[,-33] 
Databis=Databis[,-30]
Databis=Databis[,-30]
Databis=Databis[,-30]
#Databis=Databis[,-1] #On enlève les appliances 

step.model.forecastf=predict(step.modelf, newdata=Databis)
step.model.forecastb=predict(step.modelb, newdata=Databis)
step.model.forecast=predict(step.model, newdata=Databis)

source("function/rmse.R")
#On peut comparer les trois rmse en les méthodes de recherche
#rmse(y=train$Appliances[level:13964],ychap=step.model.forecastf) : 80.8521408523872
#rmse(y=train$Appliances[level:13964],ychap=step.model.forecastb) : 80.8443308079354
#rmse(y=train$Appliances[level:13964],ychap=step.model.forecast) : 80.8443308079354

plot(train$Appliances[level:13964], type='l',ylab='original',xlab='temps')
lines(step.model$fitted.values, col='red')
Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select

On cherche désormais à ajuster le nombre de degrés de libertés du modèle de la régression linéaire. On se propose tout d'abord un petit code pour voir le nombre de degrés de liberté minimisant le rmse sur un échantillon proche des valeurs à prédire.

In [24]:
rmsef=rep(0,40)#creation d'un vecteur initial de 
for (i in 0:39)
{
  step.modelf=stepAIC(fitstart, direction = "forward",scope=formula(fitall),trace=FALSE,steps=i)
  step.model.forecastf=predict(step.modelf, newdata=Databis)
  rmsef[i+1]=rmse(y=train$Appliances[level:13964],ychap=step.model.forecastf)
}
c=seq(0,40)
plot(rmsef,type='l')

On cherche donc à calculer la trace de la matrice X(XX)1X. On obtient alors 12 degrés de libertés, l'on choisit donc de prendre 12 varaiables explicatives.

In [62]:
source("function/rmse.R")
X=as.matrix(Data[-1])
M=X%*%(ginv(t(X)%*%X))%*%t(X) #On cherche à connaitre le nombre de degré de liberté 
deg=sum(diag(M))#le nombre de degrés de libertés est 12

step.modelf=stepAIC(fitstart, direction = "forward",scope=formula(fitall),trace=F,steps=deg)
step.model.forecastf=predict(step.modelf, newdata=Databis)
#rmse(y=train$Appliances[level:13964],ychap=step.model.forecastf) le rmse est de 78.9374621296323

plot(train$Appliances[level:13964],type='l')
lines(step.model.forecast,col='red')

On remarque que les données ont du mal à absorber les pics. La regression linéaire à du mal à prédire les valeurs extrêmes.

Remarque : On voit clairement que notre prédiction peut être améliorée

In [63]:
fitstart=lm(Appliances~1,data=Data)
fitall=lm(Appliances~.,data=Data)

step.model <- stepAIC(fitstart,scope=formula(fitall), direction = "forward",data=Data,steps=12,trace=F)


step.model.forecast <- predict(step.model, newdata=test[(first+1):5771,])
res[(first+1):5771]=step.model.forecast

ggplot()+geom_line(aes(as.numeric(test$date),res))

Finalement on a un score Kaggle - RMSE : 67.9

Amélioration du modèle linéaire par classification

Une fois notre meilleur modèle linéaire choisit on décide de l'améliorer en classifiant les données. Nous montrons plus tard l'importance de certaines données, notamment la température de la cuisine et son humidité (qui indiquent clairement que la famille est en train de préparer le repas par exemple) où la température et l'humidité de la salle de bain. Si ces varaiables ont des valeurs fortes nous sommes sur des pics pour appliances et ainsi nous boostons les pics de notre régression linéaire à l'aide d'une méthode proche de celle des K plus proche voisins.

Notre code comporte 4 étapes :

  • si les habitants dorment, le modèle prend des valeurs de consommation faibles;
  • si les habitants cuisinent, on choisit de placer un pic;
  • si les habitants se lavent, on choisit de placer un pic;
  • si il ya du repassagen on place un pic;
  • si les habitants partent en week end, le modèle prend des valeurs de consommation plus faibles.
In [68]:
data=total[15000:18000,]

#La cuisine est elle utilisée ?
ind=which(data$T1>25.6)  
ind2=which(test$T1[4655:5771]>25.55)  
pic4=mean(data$Appliances[ind])
res[4655:5771][ind2][5]=pic4

#La douche ou baignoire est elle utilisée ?
ind=which(data$RH_5>70 & data$T5>25.1)          
ind2=which(test$RH_5[4655:5771]>70 & test$T5[4655:5771]>25)          
pic5=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic5

#Est ce qu'un habitant fait du repassage ?
ind=which(data$RH_7>40 & data$T7>25.1 & data$NSM>=63000 & data$NSM<=63600)          
ind2=which(test$RH_7[4655:5771]>40 & test$T7[4655:5771]>24.5)          
pic7=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic7

### Si il fait beau le weekend ils vont dehors
ind=which(data$lights==0 & data$WeekStatus=='Weekend' & data$T_out>23) 
ind2=which(test[4655:5771,]$lights==0 & test[4655:5771,]$WeekStatus=='Weekend' & test[4655:5771,]$T_out>23)          
picw=mean(data$Appliances[ind])
res[4655:5771][ind2]=picw

#Les gens dorment ils ?
ind=which(data$NSM>84600) 
ind2=which(0<test[4655:5771,]$NSM & test[4655:5771,]$NSM<3600)          
pic1=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic1

ind=which(0<data$NSM & data$NSM<3600) 
ind2=which(0<test[4655:5771,]$NSM & test[4655:5771,]$NSM<3600)          
pic1=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic1

ind=which(3600<=data$NSM & data$NSM<14000) 
ind2=which(3600<=test$NSM[4655:5771] & test$NSM[4655:5771]<14000)          
pic2=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic2

ind=which(14000<=data$NSM & data$NSM<24600) 
ind2=which(14000<=test$NSM[4655:5771] & test$NSM[4655:5771]<24600)          
pic3=mean(data$Appliances[ind])
res[4655:5771][ind2]=pic3

ggplot()+geom_point(aes(x=as.numeric(test$date[4655:5771]),y=step.model.forecast),size=0.1)+geom_line(aes(as.numeric(x=test$date[4655:5771]),y=res[4655:5771]))+xlab('date')+ylab('Prédiction')

Finalement on a un score Kaggle - RMSE : 67.8

Modèle GAM

Choix du modèle GAM

Nous avons choisi de faire un modèle GAM. Nous souhaitons savoir quelle variable explicative mettre dans le modèle et confimer nos intuitions.

In [75]:
library(corrplot)
library(ellipse)
library(qgraph)
corr=data.frame(Appliances=train$Appliances,NSM=train$NSM,lights=train$lights,T1=train$T1,RH_1=train$RH_1,RH_5=train$RH_5,T5=train$T5)
corr=as.matrix(corr)
qgraph(cor(corr))

On décide alors de choisir en partie NSM et lights comme variables explicatives (bien que deux variables corrélés n'implique que l'une explique l'autre) pour commencer. Nous avons ensuite complété notre choix de variables avec un StepGam. Voici en noir la courbe reélle d'une partie des valeurs de test ayant servi d'entrainement et en rouge la prédiction de ses valeurs par le gam.

In [71]:
### Step.Gam ###
library(gam)

Gam.object=gam(step.modelf,data=train)
step.object <-step.Gam(Gam.object,trace=F,direction='forward',scope=list("lights"=~1+lights+s(lights),"RH_1"=~1+RH_1+s(RH_1),"RH_2"=~1+RH_2+s(RH_2),"RH_3"=~1+s(RH_3),"RH_4"=~1+s(RH_4),"RH_5"=~1+s(RH_5),"RH_6"=~1+s(RH_6),"RH_7"=~1+s(RH_7),"RH_8"=~1+s(RH_8),"RH_9"=~1+s(RH_9),"RH_out"=~1+RH_out+s(RH_out),"Visibility"=~1+Visibility+s(Visibility),"Posan"=~1+Posan+s(Posan),"NSM"=~1+NSM+s(NSM),"T1"=~1+T1+s(T1),"T2"=~1+T2+s(T2),"T6"=~1+T6+s(T6),"T7"=~1+T7+s(T7),"T_out"=~1+T_out+s(T_out),"BE_load_actual_entsoe_transparency"=~1+BE_load_actual_entsoe_transparency+s(BE_load_actual_entsoe_transparency),"BE_load_forecast_entsoe_transparency"=~1+BE_load_forecast_entsoe_transparency+s(BE_load_forecast_entsoe_transparency)))
Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"
In [73]:
source("function/rmse.R")
predict_gam=predict(step.object,train[10000:13964,])
# rmse(train$Appliances[10000:13964],predict_gam) de 84
plot(train$Appliances[10000:13964],type='l')
lines(predict_gam,col='red')

Finalement on a un score Kaggle - RMSE : 70

Modèle mixte : GAM et Linéaire

Nous avons améliorer nos deux modèles de prédiction GAM et linéaire en faisant une combinaison de ces mmodèles. L'idée est de laisser le modèle GAM s'exprimer lorsqu'il y a des pics importants et lorsque les habitants dorment (pour combler les défaults).

On effectue donc tout d'abord une prédiction (gt) avec un modèle gam expliquant les Appliances uniquement par NSM. On souhaite avoir un modèle prédictif valant zéro lorsque les habitants dorment et valant 1 lorsqu'il y a un pic de consommation. On pose donc Gt=gtmin(gt)max(gt)min(gt).

Nos Appliances prédites seront alors (At) telque : At=51(1Gt)+LtGt

Lt correspond à notre meilleur prédiction linéaire des Appliances manquantes du test.

Justification de ce modèle : Lorsque les habitants dorment, Gt est proche de 0 donc At équivaut à 51 ce qui est en moyenne la valeur des Appliances de minuit à 7h00 du matin. Quand il y a un pic Gt est proche de 1 et donc At correspon à la prédiction linéaire Lt.

In [89]:
#On nettoye les données pour la régression linéaire
Data=train[,-1]
Data=Data[,-33]
Data=Data[,-30]
Data=Data[,-30]
Data=Data[,-30]

#On effectue la regression linéaire
fitstart=lm(Appliances~1,data=Data)
fitall=lm(Appliances~.,data=Data)

step.model <- stepAIC(fitstart,scope=formula(fitall), direction = "forward",data=Data,steps=12,trace=F)

step.model.forecast <- predict(step.model, newdata=test[4655:5771,])
res[4655:5771]=step.model.forecast
model=gam(Appliances~s(NSM),data=train)
p=predict(model,newdata=test[4655:5771,])
p=(p-min(p))/(max(p)-min(p))

#On effectue le modèle mixte
res2=rep(1,(5771-(4655)))
res2=51.2*(1-p)+res[(first+1):5771]*p
mix=data.frame(Appliances=res2)

ggplot()+geom_point(aes(as.numeric(test$date[4655:5771]),step.model.forecast),size=0.1)+geom_line(aes(as.numeric(test$date[4655:5771]),mix$Appliances))+xlab('date')+ylab('Prédiction')
Warning message in model.matrix.default(mt, mf, contrasts):
"non-list contrasts argument ignored"

Finalement on a un score Kaggle - RMSE : 67.9

Utilisation d'un GBM

On peut visualiser un indicateurs de la performance : l'erreur Out Of Bag qui doit normalement converger vers 0 quand on augmente le nombre d'itération dans l'algorithme GBM.

Definition : L'Erreur Of Bag est E((m(X)Y)2) ou m(X) est la prédiction en fonction des variables explicatives X et Y la variable cible (ici Appliances). L'Erreur OOB représente l'écart moyen entre la prédiction et la variable cible, en régression on choisit naturellement pour formule 1nk=1n(Y^kYk)2 avec Y^=m(X).

#On réduit le jeu de données, on garde celle qui nous interresse et on enleve les redondances :
library(dplyr)
train_gbm=(train %>% select(-c(date,Day_of_week,DayType,lights,Visibility,InstantF,Posan,Month,Heure,BE_load_actual_entsoe_transparency,
                               BE_load_forecast_entsoe_transparency,BE_wind_onshore_capacity,BE_wind_onshore_capacity,BE_wind_onshore_generation_actual,BE_wind_onshore_profile,WeekStatus)))



library(gbm)
object_gbm=gbm(formula =my_formula ,

               data =train_gbm, distribution = "gaussian",
               , var.monotone = NULL, n.trees = 100,
               interaction.depth = 1, n.minobsinnode = 10, shrinkage = 0.05,
               bag.fraction = 0.5, train.fraction = 0.5, cv.folds = 20,
               keep.data = TRUE, verbose = TRUE,   n.cores = 4)

1.png

On remarque bien que certaines variables ont une influence relative faible. En appliquant le principe de parcimonie, nous avons donc décidé de retirer ces variables pour ne garder que les plus significatives.

2.png

Nous remarquons une similitude : ce sont les mêmes variables qui sont retenues dans le cadre de GBM et dans le cadre du modèle linéaire obtenu par step.AIC. On peut donc entrainer un modèle linéaire obtenu par step.AIC à l'aide d'un GBM.

Utilisation de l'interface XGBoost (Xtreme Gradient Boosting)

Nous avons utilisé XGBoost en tant que GBM amélioré.En effet,XGBoost tout comme GBM utilisent le principe de gradient Boosting, mais il nous a semblé que XGBoost proposait une implémentation plus performante que GBM, grâce à notamment:

  • parallélisation du calcul

  • tous les paramètres clés de l'apprentissages sont modifiables

  • utilisation de matrices pleines pour accélerer le temps de calcul.

  • La modélisation a été concue de telle sorte à limiter l'overfitting

En fait, nous avons utilisé XGBoost comme une interface, qui regroupe en fait plusieurs modèles. Premièrement, on définit une structure de contrôle qui permet de choisir la méthode de validation : nous avons successivement choisi la validation croisée (cv) et l'Erreur Out of Bag (OOB).

Nous nous sommes alors demandé comment les arbres de régression développés dans la phase d'apprentissage allaient mesurer l'homogénéité dans les données de Appliances, variable quantitative.

##importation des différentes bibliothèques
library(tidyverse)
library(lubridate)
library(xts)
library(dygraphs)
library(ranger)
library(dplyr)

library(mltools)
library(data.table)
library(xgboost)
library(caret)

p2=read.csv(file="submission_73.2.csv", sep=",", dec='.')$Appliances[4655:5771]
p3=read.csv(file="submission_73.3.csv", sep=",", dec='.')$Appliances[4655:5771]
p4=read.csv(file="submission_73.9.csv", sep=",", dec='.')$Appliances[4655:5771]

p_moy=(p2+p3+p4)/3
p_moy_original=p_moy

p_moy=list(p_moy)
###Méthode de complétion partielle pour conserver les valeurs d'Appliances déjà prédites et vraisemblables
p_moy_df=setDT(p_moy, keep.rownames=FALSE, key=NULL, check.names=FALSE)
p_moy$to_use=rep(0, times = 1117)
p_moy$to_use[which(p_moy$V1<80)]=1
p_moy$V1[which(p_moy$to_use==0)]=NA
par(mfrow=c(2,2))

plot(ts(p_moy_original))
lines(ts(p_moy$V1),
      col="red")

###On cree une copie de total et on l'appelle heavy pour insister sur le fait qu'il contient toutes les variables explicatives.
total_heavy=total
total_heavy_matrix=data.matrix(total_heavy, rownames.force = NA)
total_heavy_numeric= (lapply(total_heavy, function(x) as.numeric(x)))
total_heavy_numeric_df=setDT(total_heavy_numeric, keep.rownames=FALSE, key=NULL, check.names=FALSE)

library(doParallel)
library(missForest)
registerDoParallel(cores=4)

###on creee une variable total_light qui contient uniquement les variables pertinentes.
total_light=(total_heavy_numeric_df %>% select(-c(rv1,rv2,lights,Tdewpoint,RH_6,Visibility,Instant,DayType,BE_load_actual_entsoe_transparency,
                                                                      BE_load_forecast_entsoe_transparency,BE_wind_onshore_capacity,BE_wind_onshore_capacity,BE_wind_onshore_generation_actual,BE_wind_onshore_profile,Month,Windspeed,Press_mm_hg,InstantF,T2,T3,RH_2,RH_8,Day_of_week,Posan)))
permutation=sample(1:19735)
total_light_permuted=total_light[permutation]

total_with_imputation=missForest(total_light, maxiter = 10 ,ntree = 10, variablewise = TRUE,
                                decreasing = FALSE, verbose = TRUE,

                                mtry =sqrt(ncol(total_light)),
                                replace = TRUE,
                                classwt = NULL, parallelize = 'forests')

##On recupere les donnees avec les Appliances imputees
data_complete=total_with_imputation$ximp
#######Approche Xgboost

library(xgboost) 

#####Définition de la structure de controle
xgb_trcontrol = trainControl(method = "cv", number = 1000,allowParallel = TRUE, verboseIter = TRUE, returnData = FALSE)

#####definitiion du grid
xgbGrid <- expand.grid(nrounds =1,  max_depth =50,colsample_bytree = 0.9, eta = 0.1,gamma=0,min_child_weight = 1,subsample = 0.5)

total_light$is_in_test[which(is.na(total_light$Appliances))]=2
to_train=total_light[which(total_light$is_in_test<2)]
y_train=total_light$Appliances[which(total_light$is_in_test<2)]
total_light_wt_App=(to_train %>% select(-Appliances))
X_train=(xgb.DMatrix(as.matrix(total_light_wt_App ),label=y_train))


xgb_model=xgb.train(params=list(booster="gbtree"),data=X_train, trControl = xgb_trcontrol, tuneGrid = xgbGrid,  method = "xgbTree",nrounds=100000)
rm(X_test)
X_test=total_light[which(total_light$is_in_test==2)]
X_test=(X_test %>% select (-Appliances))
X_test=xgb.DMatrix(as.matrix(X_test))
predictions=predict(xgb_model,newdata=X_test)
plot(ts(predictions))

reconstitution_test=test_total$Appliances
length(which(is.na(reconstitution_test)))
tokaggle=c(appliances_imputed,predictions)
plot(tokaggle)

total_light$Appliances[which(is.na(total_light$Appliances))]=predictions
length(which(is.na(test_total$Appliances)))
test_total$Appliances[which(is.na(test_total$Appliances))]=predictions
which(is.na(test_total$Appliances))
plot(test_total$Appliances[4655:5771])
which(is.na(total_light$Appliances))

tokaggle=test_total$Appliances[4655:5771]
plot(to_g)

Techniques de prévision hybride, imputation des valeurs manquantes vs prevision pure

On remarque que le test a été echantilloné de manière linéaire dans le train. On dispose donc d'un échantillonage de bonne qualité, car les modèles peuvent être entrainés sur le train, ce qui aurait été plus difficile si les valeurs manquantes de Appliances avaient été réparties sans homogénéité.

3.png

Voici le résultat de l'imputation des données manquantes à l'aide d'une méthode de prévision à l'aide de Xgboost en utilisant une méthode de validation croisée.

4.png

Par ailleurs, voici le résultat de l'imputation des données manquantes à l'aide d'une méthode de complétion par forêts aléatoires. Plus précisement, nous avons utilisé le packagemissForest. Dans cette méthode d'imputation, l'erreur Out of Bag converge vers 0, ce qui indique le validité de notre imputation.

5.png

Finalement on a un score Kaggle - RMSE : 69.1

Remarque : Nous choisissons de travailler sur les jeu de données entier : on considère que les données imputées font partie des données initiales, donc on peut les réutiliser pour entrainer nos modèles. Nous allons discuter de ce choix dans la suite de notre étude.

Tests de nos modèles à l'aide des données imputées

Précedemment, nous avons établi que les données imputées étaient très proches des vraies données au sens du RMSE.Dès lors, on fait l'hypothèse que ces données sont plausibles.On peut donc considérer ces données manquantes avec deux points de vue:

1) Ces données manquantes peuvent être complétées par imputation (ACP par missMDA, forêts aléatoires par missForest)

2) Ces données peuvent être prédites à l'aide d'un modèle qui a été entrainé sans prendre en compte ces données.

6.png

Sur cette figure, on a représenté 3 séries mises en jeu dans le test

  • en vert, les valeurs imputées par missForest,

  • en violet, la partie présente dans le train,

  • en bleu,la partie de prédiction pure.

Conclusion : Pour améliorer notre score, nous décidons d'améliorer notre prédiction pure dans notre algorithme.

Nous représentons ici plusieurs prédictions, en utilisant plusieurs jeux de paramètres dans l'algorithme XGBoost

7.png

De plus, puisque tous les prédicteurs sont numériques,on se ramène en fait à une régression constante par morceaux sur des pavés de Rp,ces pavés étant déterminés par dichotomie succesives lors de la création de l'arbre de régression

Méthode de permutation-imputation

Nous avons constaté que, une fois le problème d'imputation résolu, le problème de prévision des valeurs dans le futur présente une difficulté supplémentaire: on doit prédire un un ensemble continu dans le temps de Appliances ,alors que dans le problème d'imputation, on pouvait compléter les valeurs manquantes grâce aux valeurs des voisins des Appliances manquantes, par exemple par la méthode des plus proches voisins.

Nous avons donc appliqué une permutation sur l'ensemble d'apprentissage afin d'homogénéiser la répartition des valeurs à prédire de facon uniforme.Ainsi, on transforme un problème de prédiction en un problème d'imputation.

8.png

En vert est représentée une prédiction réalisée à l'aide d'une régression linéaire classique.

En bleu est représentée une prédiction réalisée à l'aide de notre méthode de permutation-imputation.

Désormais, on a trois possibilités de prédire les valeurs manquantes:

1) Prévision à l'aide d'un modèle (régression, GBM, Xgboost, Random Forest)

2) Imputation, sans permutation.

3) Imputation, avec permutation.

9.png

Interprétation : Sur ce graphe, on superpose les résultats issus d'une même méthode d'imputation, avec et sans permutation. Avec permutation, les valeurs à prédire sont dispersées au sein du jeu de données et le lien temporel entre chaque relevés successif à predire est moindre. Ainsi la variance est plus élevée.

11.png

L'ACF de la série des Appliances d'entrainement montre bien une corrélation de la série à son passé.

Méthode d'amélioration des résultats : Complétion partielle

On dispose désormais de prévisions de bonne qualité, c'est-à-dire avec un faible RMSE. On choisit nos 3 meilleures prédictions au sens du RMSE, puis on crée la série qui correspond à la moyenne de ces prévisions. Le signal résultant est intéressant car on retrouve de faible appliances la nuit. Ainsi les valeurs d'Appliances très faibles (seuil Appliances <80) paraissent vraisemblables. Nous choisissons donc de supposer ces valeurs comme vraies. On réduit ainsi le cardinal de l'ensemble des données à prédire.En fixant ce seuil , on ne prédit par exemple que 893 valeurs au lieu de 1117. En rouge sont représentés les valeurs de Appliances qu'on introduit dans le jeu d'entrainement, et qu'on a donc plus besoin de prédire.

12.png

Conclusion

On peut comparer les avantages et les inconvenients de nos diverses méthodes. Dans un premier temps, nous avons mis en place une simple regression linéaire.Cette méthode a le mérite d'être facilement interprétable, et elle constitue une méthode qui conserve une continuité.En effet, de petites variations sur les données conduiront à de petites variations sur les résultats, grâce à la continuité du produit matriciel. Cependant, les régressions linéaires sont par nature inadaptées pour représenter des variations brusques (à savoir les pics).

Au contraire, les méthodes d'arbre permettent de modéliser de facon plus fidèle les discontinuités.Cependant, de faibles variation sur les données peuvent conduire à de grandes variations sur les résultats.Cependant, ces défauts ont été écartés grâce à l'utilisation du bootstrap.On a aussi établi que notre cas, l'utilisation d'un arbre de régression s'identifiait à une régréssion linéaire sur des fonctions constantes par morceaux, sur des pavés de Rp.

Remarque générale : Lors de l'utilisation de nos modèles, nous avons été confrontés au compromis biais-variance. Des modèles simples, avec peu de variables explicatives ont eu généralement un fort biais mais une variance peu élevée, et les modèles complets ont parfois réussi à obtenir un biais inférieur mais une variance élevée.

Nous sommes deuxième de la compétition.